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ABSTRACT:  A  three-dimensional  finite-volume  Eulerian-Lagrangian  Localized  Adjoint  Method  (ELLAM) 
has  been  developed,  tested,  and  successfully  implemented  as  part  of  the  U.S.  Geological  Survey  (USGS) 
MODLLOW/MOC3D  ground-water  modeling  package.  The  USGS  ELLAM  code  simulates  solute  transport 
in  ground  water  for  a  single  dissolved  constituent  subject  to  advective  transport,  hydrodynamic  dispersion  (in¬ 
cluding  mechanical  dispersion  and  diffusion),  mixing  from  fluid  sources,  and  simple  chemical  reactions  (in¬ 
cluding  linear  sorption  and  decay).  The  implementation  conserves  mass  locally  and  globally.  This  ELLAM  al¬ 
gorithm  incorporates  an  implicit-in-time,  finite-difference  approximation  for  the  dispersive  and  source/sink 
terms,  allowing  large  transport  time  increments  to  be  used  for  greater  efficiency.  It  uses  a  forward  tracking  ^- 
proach  to  move  mass  to  the  new  time  level,  distributing  mass  among  destination  cells  using  approximate 
characteristic  functions.  Numerous  test  cases  indicate  that  the  method  can  usually  yield  excellent  results,  even 
if  relatively  few  transport  time  steps  are  used,  although  the  quality  of  the  results  is  problem  dependent. 


1  INTRODUCTION 

The  modular  finite-difference  ground-water  flow 
model  (MODLLOW)  developed  by  the  U.S.  Geo¬ 
logical  Survey  (USGS)  is  a  widely  used  and  flexible 
computer  program  for  simulating  three-dimensional 
ground-water  systems  (McDonald  &  Harbaugh, 
1988,  Harbaugh  &  McDonald,  1996).  MOC3D  is  a 
solute-transport  program  that  is  integrated  with 
MODLLOW  and  has  the  capability  to  calculate 
changes  in  concentration  of  a  single  solute  subject  to 
advection,  dispersion,  diffusion,  fluid  sources,  de¬ 
cay,  and  retardation  (Konikow  et  al.  1996,  Kipp  et 
al.  1998).  MOC3D  solves  the  solute-transport  equa¬ 
tion  in  three  dimensions  using  the  method  of  char¬ 
acteristics,  with  forward  particle  tracking  to  repre¬ 
sent  advection,  coupled  with  either  an  explicit  or 
implicit  finite-difference  method  to  calculate  disper¬ 
sive  flux.  This  approach  is  optimal  for  advection- 
dominated  systems,  which  are  typical  of  many  field 
problems  involving  ground-water  contamination,  as 
it  minimizes  numerical  dispersion.  The  model  as¬ 
sumes  that  fluid  properties  are  homogeneous  and  in¬ 
dependent  of  concentration.  The  solution  techniques, 
however,  do  not  guarantee  a  mass  balance  and  also 
require  the  use  of  an  areally  uniform  grid. 

A  Linite- Volume  Eulerian-Lagrangian  Localized 
Adjoint  Method  (LVELLAM)  (Healy  &  Russell, 
1993)  was  developed  as  an  alternative  numerical 


solution  algorithm  for  the  MOC3D  transport  model. 
ELLAM  (Celia  et  al.  1990)  solves  a  mass- 
conservative  integral  form  of  the  solute-transport 
equation.  The  ELLAM  algorithm  uses  an  implicit 
time  method  for  dispersion  calculations,  which  al¬ 
lows  for  large  time  steps  without  stability  con¬ 
straints.  ELLAM  uses  an  Eulerian-Lagrangian  ad¬ 
vection  approach,  tracking  mass  through  time  and 
then  solving  a  dispersion  equation  on  a  fixed-in- 
space  grid.  Lor  advection-dominated  problems,  this 
has  the  advantage  of  generating  less  numerical  dis¬ 
persion  than  standard  Eulerian  approaches  using  fi¬ 
nite-difference  and  finite-element  methods.  ELLAM 
solves  integral  equations  and  thus  tracks  mass  asso¬ 
ciated  with  fluid  volumes,  so  that  it  conserves  mass 
locally  and  globally. 


2  GOVERNING  EQUATIONS 

The  ground-water  flow,  interstitial  velocity,  and  sol¬ 
ute-transport  equations  used  in  MOC3D  are  given  by 
Konikow  et  al.  (1996).  Solution  to  the  flow  equation 
provides  the  interstitial  velocity  field,  which  couples 
the  solute-transport  equation  to  the  ground-water 
flow  equation. 

The  governing  equation  for  this  finite- volume  ap¬ 
proach  is  an  integral  form  of  the  solute-transport 
equation,  which  is  a  statement  of  conservation  of 
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mass  over  the  domain  of  integration.  Integration 
against  a  “test  function”  (m)  provides  the  formulation 
of  conservation  of  mass,  including  treatment  of  cell 
or  subdomain  boundary  conditions  and  solute  decay. 
The  test  function  effectively  specifies  the  domain  of 
integration  for  the  transport  equation  by  the  portion 
of  the  space-time  domain  where  its  value  is  nonzero. 

If  we  multiply  the  solute-transport  equation  by 
the  test  function  u  and  integrate  over  time  and  space, 
we  have: 


f  f 
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where  Q.  is  the  entire  spatial  transport  subdomain,  T 
is  the  end  of  the  simulation  time  period  starting  at 
time  zero,  £  is  the  effective  porosity  (dimensionless), 
C  is  volumetric  concentration  (mass  of  solute  per 
unit  volume  of  fluid,  ML'^),  t  is  time  (T),  is  the 
retardation  factor,  V  is  the  interstitial  fluid  velocity 
(LT-'),  D  is  a  second-rank  tensor  of  dispersion  coef¬ 
ficients  (L^T"'),  IT  is  a  volumetric  fluid  sink  or 
source  rate  per  unit  volume  of  aquifer  (T^),  C'  is  the 
volumetric  concentration  in  the  sink/source  fluid 
(ML'^)  (if  IT  represents  a  sink,  C'  =  C),  and  A  is  the 
decay  rate  (T'^). 

The  Eulerian-Lagrangian  aspects  of  the  method 
derive  from  the  requirement  that  the  test  function 
satisfy  the  adjoint  equation. 
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Thus,  for  the  time  step  from  F  to  F+i, 

u{x,t)  =  f{x,t)e  ^  , 

where 
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so  that  /  is  constant  along  characteristics  of  the  re¬ 
tarded  interstitial  velocity  field.  With 


(5) 


(that  is,/=  1)  we  arrive  at  a  statement  of  global  con¬ 
servation  of  mass  for  a  time  step.  For  a  local  conser¬ 
vation  equation  on  each  finite-difference  cell  in 
the  transport  subdomain,  let 

Ui(x,t)  =  fi(x,t)e  (6) 


where  f/(x,t'^+^)  =  1  on  Q./  and  //(v,F+i)  =  0  else¬ 
where: 
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where  d  ■  signifies  the  spatial  boundary  of  the  argu¬ 
ment;  supp-  denotes  the  support  of  a  function  (that  is, 
the  part  of  its  domain  where  a  function  assumes  a 
non-zero  value);  n  is  the  unit  outward  normal  vector 
on  the  specified  boundary;  Qf  is  the  part  of  the  spa¬ 
tial  domain  holding  mass  at  time  for  destination 
cell  Q/  at  time  under  advection;  r"+i  = 

df2x(t",F+i)  is  the  space-time  boundary  at  time  step 
n-i-1;  and  dx  and  ds  signify  differential  volume  and 
area,  respectively. 

Note  that  equation  7  appears  as  space-time  inte¬ 
grals  of  the  dispersion  equations.  EELAM  can  be 
viewed  as  a  method  of  characteristics,  tracking  mass 
along  streamlines  of  the  flow  to  accumulate  data  to 
the  right-hand  side  of  the  system  of  equations. 


3  NUMERICAL  METHODS 
3 . 1  ELLAM  overview 

A  numerical  solution  of  the  three-dimensional 
ground-water  flow  equation  is  obtained  by  the 
MODFLOW  code  using  implicit  (backward-in-time) 
finite-difference  methods.  After  the  head  distribution 
has  been  calculated  for  a  given  time  step  or  steady- 
state  flow  condition,  the  specific  discharge  across 
every  face  of  each  finite-difference  cell  within  the 
transport  subdomain  is  calculated  using  a  finite- 
difference  approximation  (see  Konikow  et  al.  1996). 
The  seepage  velocity  is  calculated  at  points  within  a 
finite-difference  cell  based  on  linearly  interpolated 
estimates  of  specific  discharge  at  those  points  di¬ 
vided  by  the  effective  porosity  of  the  cell  (see 
Konikow  et  al.  1996). 

Advection  in  flowing  ground  water  is  simulated 
by  mass  tracking  along  the  characteristic  curves  de¬ 
termined  by  the  seepage  velocity.  Calculation  of  ad- 
vective  movement  during  a  flow  time  step  is  based 
on  the  specific  discharges  computed  at  the  end  of  the 
time  step. 

As  in  MOC3D,  tracking  is  performed  using  lin¬ 
ear  interpolation  of  velocity  in  the  direction  of  the 
component  of  interest  and  piecewise-constant  inter- 


polation  in  the  other  two  directions.  The  approach  is 
to  solve  a  system  of  three  ordinary  differential  equa¬ 
tions  to  find  the  characteristic  curves  [x  =  x(t),  y  = 
y(t),  and  z  =  z(tj\  along  which  fluid  is  advected.  This 
is  accomplished  by  introducing  a  set  of  moving 
points  that  can  be  traced  within  the  stationary  coor¬ 
dinates  of  a  finite-difference  grid.  Each  point  corre¬ 
sponds  to  one  characteristic  curve,  and  values  of  v,  y, 
and  z  are  obtained  as  functions  of  t  for  each  charac¬ 
teristic  (Garder  et  al.  1964).  Each  point  moves 
through  the  flow  field  at  a  rate  governed  by  the  flow 
velocity  acting  along  its  trajectory. 

The  EEEAM  equations  suggest  that  mass  is 
tracked  backwards  along  characteristics  to  the  pre¬ 
image  of  each  cell  or  boundary  face.  It  is  not  possi¬ 
ble,  however,  to  exactly  locate  all  of  the  mass  at  the 
previous  time  level  by  backtracking  a  finite  number 
of  points.  In  order  to  achieve  mass  balance,  this  im¬ 
plementation  of  the  EEEAM  algorithm  tracks  the 
known  mass  distribution  forward  from  the  old  time 
level  to  the  new  time  level. 

Numerically,  this  is  accomplished  by  tracking 
mass  forward  from  the  old  time  level,  n,  along  char¬ 
acteristics.  Each  cell  is  divided  into  subcells,  where 
the  number  of  subcells  is  determined  by  parameters 
NSC,  NSR,  NSE,  specifying  the  number  of  subcells 
in  the  column,  row,  and  layer  direction,  respectively. 
The  center  of  each  subcell  is  the  location  at  time 
level  n  of  one  of  the  moving  points  discussed  above, 
which  is  tracked  through  the  time  step  under  advec- 
tion.  Depending  on  the  exact  location  of  this  point  in 
the  destination  cell  at  the  new  time,  all  of  the  mass 
in  the  subcell  may  or  may  not  also  be  found  in  that 
destination  cell.  In  order  to  mitigate  the  effects  of 
unwarranted  mass  lumping,  subcell  mass  is  distrib¬ 
uted  among  cells  neighboring  the  destination  cell 
using  “approximate  test  functions,”  W;.  The  value  of 
Wi  at  the  subcell  center  destination  point  is  the  frac¬ 
tion  of  subcell  mass  to  be  distributed  to  cell  /. 

This  yields  the  formulation  of  the  time  step  n 
term  in  equation  7, 
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where  summation  runs  through  all  subcells  of  each 
cell  in  the  transport  subdomain,  b  is  the  cell  thick¬ 
ness,  and  pf  is  the  image  of  p  under  forward  tracking 
to  the  new  time  level. 

The  model  assumes  a  source  or  sink  acts  uni¬ 
formly  over  the  entire  cell  surrounding  a  source  or 
sink  node.  Eor  a  cell  containing  a  fluid  source,  a  sin¬ 
gle  time  step  is  discretized  into  a  number  of  sub-time 
steps  determined  by  parameter  NT,  and  the  com¬ 
posite  trapezoidal  integration  rule  is  applied  in  time. 
Mass  is  tracked  to  varying  locations  within  the 


transport  subgrid  depending  on  when  in  the  transport 
time  step  the  mass  enters  the  system.  At  each  sub¬ 
time  step,  inflow  mass  is  spatially  discretized, 
tracked,  and  accumulated  just  like  mass  that  began 
the  time  step  already  in  the  system,  but  for  the 
shorter  interval. 

Sink  concentration  is  assumed  to  be  the  average 
nodal  concentration  for  the  transport  time  step,  with 
the  exception  of  a  sink  due  to  evapotranspiration, 
where  sink  concentration  is  taken  to  be  zero.  Inte¬ 
gration  rules  are  midpoint  in  space  and  a  one  point 
backward  Euler  in  time. 

The  quantity  mass/porosity  in  a  cell  at  the  new 
time  level  F+i  is  expressed  using  the  trapezoidal  rule 
for  integration,  formulated  over  each  cell  octant. 
Concentrations  at  octant  corners  are  weighted  aver¬ 
ages  of  neighboring  node  concentrations,  determined 
by  trilinear  interpolation. 

Eike  the  algorithms  in  the  previous  versions  of 
MOC3D,  the  EEEAM  method  approximates  total 
solute  flux  across  the  transport  subdomain  boundary 
by  the  advective  flux.  This  approximation  is  not  de¬ 
manded  by  EEEAM  methods  in  general,  but  it  is  a 
feature  of  this  particular  implementation.  This  ap¬ 
proximation  means  that  boundary-face  concentra¬ 
tions  are  not  coupled  to  cell-center  concentrations 
through  the  numerical  derivative.  All  mass  moving 
in  and  out  of  the  transport  subdomain  can  be  tracked 
using  the  advective  algorithm.  Mass  tracked  across 
outflow  boundaries  provides  data  for  a  system  of 
outflow  boundary  equations  decoupled  from  the  cell 
equations.  User-input  inflow  concentrations  together 
with  the  outflow  solutions  then  accumulate  on  the 
right-hand  side  of  the  system  of  cell  equations. 

Eor  an  inflow  boundary,  as  for  a  source,  a  single 
time  step  is  discretized  into  a  number  of  sub-time 
steps  determined  by  parameter  NT.  The  only  differ¬ 
ence  in  the  treatment  of  the  inflow  boundary  from 
the  treatment  of  the  source  is  that  only  the  two- 
dimensional  boundary  face  is  discretized,  while  for  a 
source,  the  entire  cell  is  discretized.  Eor  a  cell 
the  integration  is  performed  over  the  intersection  of 
the  support  of  the  space-time  test  function  for  that 
cell  and  the  transport  subdomain  boundary;  that  is, 
all  mass  entering  through  the  boundary  and  advected 
to  Q-i  during  the  time  step  is  accumulated  to  the 
right-hand  side  of  local  equation  1. 

The  right-hand-side  outflow  boundary  integrals 
are  constructed  from  the  mass  contributions  tracked 
across  the  boundary  from  interior  cells,  sources,  and 
inflow  boundaries  during  the  transport  time  step.  All 
mass  associated  with  a  tracked  point  that  reaches  the 
outflow  boundary  at  any  time  during  the  time  step  is 
considered  to  leave  the  transport  subdomain.  Test 
functions  are  evaluated  to  distribute  mass  among 
neighboring  boundary  outflow  faces. 


3.2  Accuracy  criteria 

An  accuracy  criterion  incorporated  in  the  model  con¬ 
strains  the  distance  that  solute  mass  is  advected 
during  each  transport  time  step.  A  restriction  can  be 
placed  by  the  model  user  on  the  size  of  the  time  step 
to  ensure  that  the  number  of  grid  cells  a  point  moves 
in  the  v-,  y-,  or  z-directions  does  not  exceed  some 
maximum.  This  translates  into  a  limitation  on  the 
transport  time-step  length.  If  the  time  step  used  to 
solve  the  flow  equation  exceeds  the  time  limit,  the 
flow  time  step  will  be  subdivided  into  an  appropriate 
number  of  equal- sized  smaller  time  increments  for 
solving  transport. 

For  advective  transport,  a  mesh  density  sufficient 
so  that  at  least  four  grid  nodes  are  represented  across 
a  solute  front  (or  zone  of  relatively  steep  concentra¬ 
tion  gradient)  is  needed  for  good  accuracy.  Similarly, 
for  advecting  a  peak  concentration,  the  area  of  the 
peak  should  be  represented  across  at  least  eight 
nodes  of  the  grid  for  good  accuracy.  In  such  cases, 
testing  suggests  that  a  peak  concentration  value  can 
be  advected  with  a  very  small  dissipation  of  the 
maximum  concentration  per  time  step  for  a  variety 
of  Courant  numbers.  With  insufficient  mesh  density, 
a  peak  will  dissipate  (or  decay)  rapidly  for  an  initial 
period  of  time  during  which  it  spreads  out  and  os¬ 
cillates;  thereafter,  the  numerical  decay  is  very  slow 
and  the  oscillations  do  not  worsen.  A  fine  discreti¬ 
zation  of  tracked  mass  (large  NSC,  NSR,  NSL)  re¬ 
duces  the  rate  of  peak  decay  when  modeling  with 
many  transport  time  steps.  Regardless  of  the  solution 
accuracy,  global  mass  is  conserved. 

The  accuracy  of  the  dispersion  calculation  is  gov¬ 
erned  in  part  by  the  accuracy  of  the  central- 
difference  approximations  to  the  space  derivatives, 
meaning  a  finer  mesh  will  result  in  better  accuracy. 
The  implicit  formulation  for  the  solution  of  the  dis¬ 
persion  equation  is  unconditionally  stable.  This  al¬ 
lows  for  large  time  steps  during  the  simulation.  Be¬ 
cause  ELLAM  solves  the  dispersion  equation  along 
characteristics,  thus  avoiding  large  values  of  the  sec¬ 
ond  time  derivative  of  the  solution  at  passage  of  a 
steep  front,  error  in  calculation  of  the  time  derivative 
may  be  expected  to  be  small  compared  to  a  standard 
finite-difference  solution  to  an  advection-diffusion 
equation.  Some  dependence  of  accuracy  of  the  dis¬ 
persion  calculation  on  the  size  of  the  time  step  is  re¬ 
tained,  however.  Note  that  stability  does  not  imply 
accuracy;  accuracy  of  the  solution  to  the  dispersion 
equation  decreases  as  the  time  step  size  increases. 
On  the  other  hand,  modeling  with  many  time  steps  in 
order  to  resolve  dispersion  to  the  desired  accuracy 
could  result  in  a  loss  of  peak  to  numerical  dispersion 
inherent  in  the  treatment  of  advection,  an  effect  that 
can  be  reduced  by  increasing  NSC,  NSR,  and  NSL. 

One  additional  difficulty  encountered  with  im¬ 
plicit  temporal  differencing  results  from  the  use  of  a 


symmetric  spatial  differencing  for  the  cross-product 
terms  of  the  dispersion  tensor.  This  creates  a  poten¬ 
tial  for  overshoot  and  undershoot  in  the  calculated 
concentration  solution,  particularly  when  the  veloc¬ 
ity  field  is  oblique  to  the  axes  of  the  grid.  A  remedy 
for  excessive  overshoot  and  undershoot  is  to  refine 
the  finite-difference  mesh.  This  may,  however,  in¬ 
crease  simulation  times. 


4  MODEL  TESTING  AND  EVALUATION 

The  ELLAM  simulator  was  tested  and  evaluated  by 
running  a  suite  of  test  cases  (see  Konikow  et  al. 
1996,  Kipp  et  al.  1998).  This  suite  includes  base  re¬ 
sults  generated  by  analytical  solutions  and  by  other 
numerical  models.  It  spans  a  range  of  conditions  and 
problem  types  so  that  the  user  will  gain  an  apprecia¬ 
tion  for  both  the  strengths  and  weaknesses  of  this 
particular  code.  All  test  cases  involve  steady-state 
flow  conditions. 

4.1  One -dimensional  flow 

The  first  test  case  evaluates  ELLAM  for  a  relatively 
simple  system  involving  one-dimensional  solute 
transport  in  a  finite-length  aquifer  having  a  third- 
type  source  boundary  condition.  The  numerical  re¬ 
sults  are  compared  to  an  analytical  solution  by 
Wexler(1992,  p.  17). 

The  length  of  the  system  is  12  cm;  other  parame¬ 
ters  are  summarized  in  Table  1.  The  solute-transport 
equation  was  solved  using  ELLAM  on  a  120-cell 
subgrid  to  assure  a  constant  velocity  within  the 
transport  domain  and  to  allow  an  accurate  match  to 
the  boundary  conditions  of  the  analytical  solution. 


Table  1.  Parameters  used  in  ELLAM  simulation 
of  solute  transport  in  a  one-dimensional,  steady- 
state  flow  system. 


Parameter  Value 


Transmissivity  {Txx) 

0.01  cm^/s 

£ 

0.1 

Longitudinal  dispersivity  (CCf) 

0.1  cm 

Total  simulated  time 

120  s 

Lx 

0.1  cm/s 

II 

0.0  cm/s 

Initial  concentration  (Cp) 

0.0 

Source  concentration  (C' ) 

1.0 

Number  of  rows  &  layers 

1 

Number  of  columns 

122 

Ax  =  Ay 

0.1  cm 

Az 

1.0  cm 

NSC 

4 

NSR=NSL 

2 

NT 

128 

Using  the  explicit  finite-difference  solution  in 
MOC3D,  the  dispersion  coefficient  imposed  the 
limiting  stability  criteria,  and  2401  time  increments 
were  required  to  solve  the  transport  equation.  The 
implicit  solver  of  MOC3D,  however,  required  only 
241  moves.  ELLAM  results  for  241  transport  time 
steps,  NSC  =  4,  NSR  =  NSL  =  2,  and  NT  =  128  are 
essentially  identical  to  the  analytical  results.  There¬ 
fore,  these  solutions  are  not  plotted.  Instead,  results 
are  plotted  for  substantially  fewer  time  increments. 

A  variety  of  parameter  values  were  evaluated. 
Figure  1  shows  the  analytical  solution  and  two  EL¬ 
LAM  solutions  for  a  low-dispersion  case.  For  121 
transport  time  increments,  using  NSC  =  32,  NSR  = 
NSL  =  2,  and  NT  =  128,  there  is  a  very  close  match 
between  the  numerical  and  analytical  solutions.  To 
improve  clarity  in  showing  the  results  for  this  case, 
only  every  fourth  data  point  is  shown,  except  for 
times  less  than  10  seconds  at  x  =  0.05,  where  every 
point  is  plotted.  The  efficiency  of  the  numerical  so¬ 
lution  can  be  improved  by  about  a  factor  of  four  by 
setting  NSC  =  4;  the  results  are  very  similar,  so  are 
not  plotted,  although  the  concentration  is  just 
slightly  low  at  the  first  grid  cell.  The  simulation  took 
261  seconds  when  NSC  =  32,  but  only  60  seconds 
when  NSC  =  4.  (Simulations  were  executed  on  a 
Data  General  Unix  workstation.)  For  12  time  incre¬ 
ments,  using  NSC  =  4,  NSR  =  NSL  =  2,  and  NT  = 
128,  concentrations  at  early  times  and  short  dis¬ 
tances  are  somewhat  low,  but  elsewhere  the  results 
look  excellent.  Thus,  there  is  an  overall  good  agree¬ 
ment  with  the  analytical  results,  as  well  as  with  the 
previously  published  MOC3D  results  that  used  20 
times  as  many  time  increments.  The  simulation  us¬ 
ing  12  time  increments  took  only  15  seconds. 


Figure  1.  Numerical  and  analytical  solutions  at  three  different 
locations  for  solute  transport  in  a  one-dimensional,  steady  flow 
field. 


In  all  cases  described  above,  the  mass-balance  er¬ 
ror  was  less  than  0.001  percent.  In  contrast,  the 
mass-balance  errors  for  these  problems  using  the  ex¬ 
plicit  and  implicit  versions  of  the  method-of- 


characteristics  code  yielded  mass-balance  errors  of 
up  to  a  few  percent  in  some  cases. 

4.2  Point  initial  condition  in  uniform  flow 

Three-dimensional  solute  transport  of  an  instantane¬ 
ous  point  source  (Dirac  initial  condition)  in  a  uni¬ 
form  flow  field  was  used  as  another  test  problem.  An 
analytical  solution  for  an  instantaneous  point  source 
in  a  homogeneous  infinite  aquifer  is  given  by  Wex- 
ler  (1992,  p.  42),  who  presents  the  POINT3  code  for 
a  related  case  of  a  continuous  point  source.  The 
POINT3  code  was  modified  to  solve  for  the  desired 
case  of  an  instantaneous  point  source. 

The  test  problem  was  designed  to  evaluate  the 
numerical  solution  for  a  case  in  which  flow  occurs  at 
45  degrees  to  the  x-  and  y-axes.  This  allows  us  to 
evaluate  the  accuracy  and  sensitivity  of  the  numeri¬ 
cal  solution  to  the  orientation  of  the  flow  relative  to 
the  grid.  The  assumptions  and  parameters  for  this 
test  case  are  summarized  in  Table  2  and  are  de¬ 
scribed  in  more  detail  by  Konikow  et  al.  (1996). 


Table  2.  Parameters  used  in  ELLAM  simulation  of 
three-dimensional  transport  from  a  point  source 
with  flow  at  45  degrees  to  x-  and  y-axes. 


Parameter 

Value 

II 

10.0  m2/day 

£ 

0.1 

CCl 

1.0  m 

i 

II 

f 

0.1  m 

Total  simulated  time 

40  days 

II 

1.0275  m/day 

0.0  m/day 

Initial  concentration  at  source 

IX  10^ 

Grid  location  of  source 

(11,36,4) 

Number  of  rows  &  columns 

72 

Number  of  layers 

24 

Ax  =  Ay 

3.33  m 

Layer  thickness  (Az) 

10.0  m 

NSC=NSR=NSL 

4 

NT 

2 

The  results  of  the  test  problem  for  flow  at  45  de¬ 
grees  to  the  grid  are  shown  in  Figure  2.  The  analyti¬ 
cal  solution  for  f  =  130  days,  which  provides  the  ba¬ 
sis  for  the  evaluation,  is  shown  in  Figure  2a.  The 
ELLAM  solution  used  the  analytical  solution  at  t  = 
90  days  as  the  initial  conditions,  so  the  elapsed  time 
for  the  comparison  is  40  days.  The  results  using 
three  transport  time  increments,  NSC  =  NSR  =  NSL 
=  4,  and  NT  =  2  are  shown  in  Figure  2b  for  the  hori¬ 
zontal  plane  of  the  initial  source.  ELLAM  produces 
the  symmetry  characteristic  of  the  analytical  solu¬ 
tion.  There  is  also  slight  longitudinal  spreading  (nu¬ 
merical  dispersion). 


The  numerical  results  in  Figure  2b  show  some 
distortion  of  the  shape  of  the  plume  relative  to  the 
analytical  solution.  It  is  not  as  pronounced,  however, 
as  the  “hourglass”  shape  yielded  by  MOC3D  for  this 
problem  (see  Kipp  et  al.  1998,  Fig.  14).  There  is  a 
narrowing  of  the  plume  calculated  with  the  numeri¬ 
cal  model,  which  is  characteristic  of  a  grid- 
orientation  effect  and  is  caused  primarily  by  the  off- 
diagonal  (cross-product)  terms  of  the  dispersion  ten¬ 
sor.  When  flow  is  oriented  parallel  to  the  grid,  or 
when  longitudinal  and  transverse  dispersivities  are 
equal,  the  cross-product  terms  of  the  dispersion  ten¬ 
sor  are  zero.  Because  flow  is  at  45  degrees  to  the 
grid  in  this  test  problem,  the  cross-product  disper¬ 
sive  flux  terms  are  of  maximum  size  and  negative 
concentrations  are  most  likely  to  occur.  The  calcu¬ 
lated  concentration  field  is  less  accurate  in  this  case 
largely  because  the  standard  differencing  scheme  for 
the  cross-product  dispersive  flux  terms  can  cause 
overshoot  and  undershoot  of  concentrations.  If  the 
base  (or  background)  is  zero  concentration,  then  un¬ 
dershoot  will  cause  negative  concentrations.  The 
magnitude  of  this  overshoot  and  undershoot  effect  is 
reduced  by  using  a  finer  grid. 


Figure  2.  Concentration  contours  for  (a)  analytical  and  (b)  EL- 
LAM  solutions  for  transport  of  a  point  initial  condition  in  uni¬ 
form  flow  at  45  degrees  to  the  x-direction  at  f  =  130  days. 
Contour  values  are  log  of  concentration. 


Indeed,  some  small  areas  of  negative  concentra¬ 
tions  were  calculated.  The  extent  of  the  areas  of 
negative  concentration  are  indicated  by  shading  all 
areas  where  the  relative  concentration  is  less  than 
-0.05  and  less  than  -10.0.  Decreasing  the  size  of  the 
transport  time  increment  did  not  substantially  reduce 
the  area  over  which  negative  concentrations  oc¬ 
curred.  The  increase  in  execution  time,  however, 
was  significant,  so  the  very  small  improvement  does 
not  appear  to  justify  the  extra  computational  costs. 

This  test  case  involving  a  Dirac  initial  condition 
with  flow  at  a  45  degree  angle  to  the  grid  represents 
a  stringent  test  for  any  solute-transport  model,  and 
the  ELLAM  results  for  this  case  are  qualitatively 
good.  Of  all  the  test  cases  for  which  ELLAM  was 
evaluated,  however,  the  results  were  least  accurate 
for  this  particular  set  of  test  conditions. 

4.3  Constant  source  in  nonuniform  flow 

Burnett  &  Erind  (1987)  used  a  numerical  model  to 
simulate  a  hypothetical  problem  having  a  constant 
source  of  solute  over  a  finite  area  at  the  surface  of  an 
aquifer  having  homogeneous  properties,  but  nonuni¬ 
form  boundary  conditions,  which  result  in  nonuni¬ 
form  flow.  Because  an  analytical  solution  is  not 
available  for  such  a  complex  system,  we  use  their  re¬ 
sults  for  this  test  case  as  a  benchmark  for  compari¬ 
son  with  the  results  of  applying  the  ELLAM  algo¬ 
rithm  in  MOC3D,  as  was  also  done  by  Konikow  et 
al.  (1996)  and  Kipp  et  al.  (1998).  Burnett  &  Erind 
(1987)  used  an  alternating-direction  Galerkin  finite- 
element  technique  to  solve  the  flow  and  solute- 
transport  equations  in  both  two  and  three  dimen¬ 
sions.  A  detailed  description  of  the  problem  geome¬ 
try  and  of  the  parameters  for  the  numerical  simula¬ 
tion  are  presented  by  Konikow  et  al.  (1996,  p.  55- 
60). 

Cases  of  both  two-  and  three-dimensional  trans¬ 
port  were  examined  for  this  problem,  but  only  the 
latter  case  will  be  presented  here.  The  grids  used  in 
the  ELLAM  simulations  were  designed  to  match  as 
closely  as  possible  the  finite-element  mesh  used  by 
Burnett  &  Erind  (1987).  Some  differences  in  discre¬ 
tization,  however,  could  not  be  avoided  because  the 
finite-element  method  uses  a  point-centered  grid 
whereas  ELLAM  uses  a  block-centered  grid.  The 
former  allows  specifications  of  values  at  nodes, 
which  can  be  placed  directly  on  boundaries  of  the 
model  domain.  Nodes  in  ELLAM  are  located  at  the 
centers  of  cells,  and  block-centered  nodes  are  always 
one-half  of  the  grid  spacing  away  from  the  edge  of 
the  model  domain.  Among  the  small  differences 
arising  from  the  alternative  discretization  schemes 
are  that,  in  the  ELLAM  grid,  (1)  the  modeled  loca¬ 
tion  of  the  14.25  m  long  source  area  is  offset  by 
0.225  m  towards  the  right,  and  (2)  the  total  length  of 
the  domain  is  199.5  m. 


The  input  data  values  for  this  analysis  are  listed  in 
Table  3.  The  top  flow  layer  consisted  of  constant- 
head  nodes  and  the  solute  source. 


Table  3.  Parameters  used  for  ELLAM  simulation 
of  three-dimensional  transport  from  a  continuous 
point  source  in  a  nonuniform,  steady-state,  flow 
system  (described  by  Burnett  &  Frind,  1987). 


Parameter 

Value 

K 

1.0  m/day 

£ 

0.35 

CCl 

3.0  m 

Uth 

0.10  m 

arv 

0.01  m 

Total  simulation  time 

12,000  days 

Source  concentration  (C' ) 

1.0 

Number  of  rows 

15 

Number  of  columns 

141 

Number  of  layers 

91 

Ax 

1.425  m 

Ay 

1.0  m 

Az 

0.2222-0.2333 

m 

NSC=NSR=NSL 

4 

NT 

32 

Results  for  the  three-dimensional  case  based  on 
the  test  case  of  Burnett  &  Frind  (1987)  are  presented 
in  Figure  3,  which  shows  the  transport  results  in  a 
vertical  plane  at  the  middle  of  the  plume.  The  EL¬ 
LAM  plume  (Lig.  3b)  closely  matches  that  calcu¬ 
lated  by  the  finite-element  model  (Lig.  3a),  although 
the  former  shows  slightly  further  downstream  mi¬ 
gration  of  low  concentrations  of  solute.  The  ELLAM 
solution  used  seven  transport  time  steps  whereas  the 
finite-element  solution  used  40  time  steps.  The  EL¬ 
LAM  solution  provides  a  closer  match  to  the  Burnett 
&  Lrind  (1987)  solution  than  do  the  previous 
MOC3D  results,  which  used  381  time  steps  with  the 
implicit  dispersion  calculation  (Kipp  et  al.  1998)  and 
4,218  with  the  explicit  dispersion  calculation 
(Konikow  et  al.  1996). 


Figure  3.  Results  of  3-D  nonuniform-flow  test  case:  (a)  finite- 
element  model  using  40  time  steps  (Burnett  &  Frind,  1987,  Fig. 
8c),  and  (b)  EFFAM  solution  using  7  time  steps.  Contours  are 
relative  concentration  (contour  interval  is  0.2). 


4.4  Model  availability 

The  new  ELLAM  code  and  a  documentation  re¬ 
port  will  be  available  for  downloading  over  the 
Internet  from  a  USGS  software  repository.  The  re¬ 
pository  is  accessible  on  the  World  Wide  Web  from 
the  USGS  Water- Resources  Information  Web  page 
at  http://water.usgs.gov/software/  or  from  an  alterna¬ 
tive  web  page  for  USGS  ground-water  models  at 
http://water.usgs.gov/nrp/gwsoftware/. 


5  RELATIVE  COMPUTATIONAL  AND 

STORAGE  ELLICIENCY 

Computer-memory  storage  requirements  for  EL¬ 
LAM  are  greater  than  those  for  the  explicit  or  im¬ 
plicit  MOC3D  dispersive  transport  algorithm.  The 
additional  arrays  required  can  increase  the  memory 
size  requirement  by  as  much  as  a  factor  of  three. 

The  computational  effort  required  by  the  ELLAM 
simulator  is  strongly  dependent  on  the  size  of  the 
problem  being  solved,  as  determined  by  the  total 
number  of  nodes,  the  NS  and  NT  values,  and  the  to¬ 
tal  number  of  time  increments.  Analyses  indicate 
that  the  greatest  computational  effort,  as  measured 
by  CPU  time,  is  typically  expended  in  the  mass 
tracking  routines.  Lor  a  given  problem,  computa¬ 
tional  time  may  vary  significantly  as  a  function  of 
the  characteristics  of  the  particular  computer  on 
which  the  simulation  is  performed,  and  on  which 
LORTRAN  compiler  and  options  were  used  to  gen¬ 
erate  the  executable  code. 


6  CONCLUSIONS 

The  accuracy  and  precision  of  the  numerical  results 
of  the  implicit  ELLAM  simulator  were  tested  and 
evaluated  by  comparison  to  analytical  and  numerical 
solutions  for  the  same  set  of  test  problems  as  docu¬ 
mented  previously  for  MOC3D,  although  the  in¬ 
stantaneous  point  source  problem  was  modified 
slightly.  These  evaluation  tests  indicate  that  the  so¬ 
lution  algorithms  in  the  ELLAM  model  can  success¬ 
fully  and  accurately  simulate  three-dimensional 
transport  and  dispersion  of  a  solute  in  flowing 
ground  water.  The  numerical  methods  used  to  solve 
the  governing  equations  have  broad  general  capabil¬ 
ity  and  flexibility  for  application  to  a  wide  range  of 
hydrogeological  problems.  To  avoid  non-physical 
oscillations  and  loss  of  peak  concentrations,  care 
must  be  taken  to  use  a  grid  having  sufficient  mesh 
density  to  adequately  resolve  sharp  fronts. 

Relative  to  the  method  of  characteristics,  the  pri¬ 
mary  advantages  of  the  ELLAM  code  are  that  fewer 
transport  time  steps  need  be  used  and  that  mass  is 
conserved  globally,  and  relative  to  other  standard 


numerical  methods,  it  minimizes  numerical  disper¬ 
sion  in  advection-dominated  systems.  Using  EL- 
LAM  with  few  time  steps  can  provide  an  accurate 
and  cost-effective  way  of  discerning  salient  features 
of  the  solute-transport  process  under  a  complex 
given  set  of  boundary  conditions.  Lurthermore,  the 
ELLAM  algorithm  eliminates  the  previous  restric¬ 
tion  in  MOC3D  that  the  transport  grid  had  to  be  uni¬ 
formly  spaced. 
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